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We investigate the phase diagram of spinless fermions with nearest and next-nearest neighbour density- 
density interactions on the honeycomb lattice at half-filling. Using Exact Diagonalization techniques of the full 
Hamiltonian and constrained subspaces, combined with a careful choice of finite-size clusters, we determine 
the different charge orderings that occur for large interactions. In this regime we find a two-sublattice Neel-like 
state, a charge modulated state with a tripling of the unit cell, a zig-zag phase and a novel charge ordered states 
with a 12 site unit cells we call Neel domain wall crystal, as well as a region of phase separation for attractive 
interactions. A sizeable region of the phase diagram is classically degenerate, but it remains unclear whether an 
order-by-disorder mechanism will lift the degeneracy. For intermediate repulsion we find evidence for a Kekule 
or plaquette bond-order wave phase. We also investigate the possibility of a spontaneous Chern insulator phase 
(dubbed topological Mott insulator), as previously put forward by several mean-field studies. Although we are 
unable to detect convincing evidence for this phase based on energy spectra and order parameters, we find an 
enhancement of current-current correlations with the expected spatial structure compared to the non-interacting 
situation. While for the studied t—\ 1 —V 2 model the phase transition to the putative topological Mott insulator 
is preempted by the phase transitions to the various ordered states, our findings might hint at the possibility 
for a topological Mott insulator in an enlarged Hamiltonian parameter space, where the competing phases are 
suppressed. 

PACS numbers: 71.10.Fd, 71.27.+a, 71.30.+h, 75.40.Mg 


I. INTRODUCTION 

The seminal discoveries of quantum Hall phases 1 have pro¬ 
vided the first examples of topological phases of matter, i.e. 
phases which cannot be adiabatically connected to conven¬ 
tional ones. More recently, several other topological insula¬ 
tors or superconductors have been theoretically predicted and 
experimentally observed 2,3 , giving rise to an intense activity 
of the community. For such materials, a very thorough un¬ 
derstanding has been achieved thanks to the fact that a non¬ 
interacting starting point is sufficient: indeed, topological in¬ 
sulators can be obtained from band theory in the presence of 
strong spin-orbit coupling. 

On the other hand, given the variety of electronic phases 
that are also generated by strong interactions, a highly topical 
question is to investigate the appearance of topological insula¬ 
tors due to correlations. This new route would not require in¬ 
trinsic spin-orbit coupling, and thus could be advantageous to 
many applications. Nevertheless, we face the well-known dif¬ 
ficulty of strongly correlated systems for which we lack unbi¬ 
ased tools, and often have to resort to some hypothesis. In this 
context, a very promising result has been obtained by Raghu 
et al. 4 who have shown the emergence of quantum anoma¬ 
lous (spin) Hall (QAH) states on the honeycomb lattice with 
strong density-density repulsions between next-nearest neigh¬ 
bors. These phases can be characterized by the existence of 
spontaneously appearing charge (respectively spin) currents 
for spinless (respectively spinful) fermions. This finding has 
also been reproduced in the spinless case using more sys¬ 
tematic mean-field approaches 5,6 . Nevertheless, in the latter 
study, 6 the QAH extension in the phase diagram has shrunk 
substantially due to the stability of charge modulation with 


larger unit cell. 

The possibility to generate topological phases with longer- 
range interactions has triggered a large theoretical activity 
to investigate its experimental feasibility using metallic sub¬ 
strate 7 , cold-atomic gases on optical lattices 8 or RKKY inter¬ 
action 9 among many others. 

In the light of several recent approximate and numerical 
studies 10-12 with conflicting conclusions regarding the pres¬ 
ence of the putative topological Mott state and competing or¬ 
dered phases, we feel that a comprehensive state-of-the-art 
exact diagonalization study could shed some light on these is¬ 
sues. Based on our simulations, we will provide numerical ev¬ 
idence that quantum fluctuations at small V 2 > 0 are responsi¬ 
ble for an enhanced response similar to the QAH state, but un¬ 
fortunately our numerical data (obtained on finite clusters up 
to N = 42 sites) indicates that these short-range fluctuations 
presumably do not turn into a true long-range ordered phase. 
In addition, we will provide a detailed analysis of the charge 
and bond orderings that occur at large interaction strength, al¬ 
lowing us to sketch a global phase diagram of the minimal 
model studied so far. In Sec. II, we detail the model and meth¬ 
ods we use. We then discuss the large interaction (classical) 
limit in Sec. III. Our numerical results are given in Sec. IV, 
and we conclude in Sec. V. 


II. MODEL, METHOD AND PHASE DIAGRAM 
A. Hamiltonian 

Following the proposal by Raghu et al. 4 of an exotic topo¬ 
logical Mott insulator in a half-filled honeycomb lattice with 
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FIG. 1. (Color online) (a) Illustration of the honeycomb lattice with 
Vi , V 2 interactions and the hopping t. (b) First (solid line) and second 
(dashed line) Brillouin zone of the honeycomb lattice including the 
location of a few special points in the Brillouin zone. 


V 2 interaction, we will consider the following Hamiltonian : 


T-L = —ty^(c\cj H- h.c.) (1) 

(ij) 

+ K i y(«,-l/2)(n r l/2) 

<*i> 

+ ( n i - 1 / 2 ) (nj - 1 / 2 ) 

«tf» 

depicted in Fig. 1(a), where c* and c\ are the spinless 
fermionic operators, t = 1 is the nearest-neighbor hopping 
amplitude, Vi and V2 are the density-density repulsion or at¬ 
traction strengths respectively on nearest- and next-nearest 
neighbors. Note that we will focus on the half-filled case 
where thanks to particle-hole symmetry, the chemical poten¬ 
tial is known to be zero exactly. 

Despite the particle-hole symmetry, it seems that quantum 
Monte-Carlo (QMC) approaches exhibit a sign problem that 
prohibits them. In fact, only in the simpler case V 2 = 0 which 
exhibits a direct phase transition from semi-metallic (SM) 
phase to charge density wave (CDW) as a function of V\ /£, 
a very interesting recent proposal was made allowing to refor¬ 
mulate the problem without sign-problem 13 , thus amenable to 
accurate, unbiased QMC simulations 14 . Quite remarkably, the 
case Vi > 0 and V2 < 0 can also be studied with QMC with¬ 
out a minus-sign problem using a Majorana representation 15 . 
In the absence of a QMC algorithm for the frustrated case, we 
use Exact Diagonalization (ED) to get unbiased numerical re¬ 
sults. Of course, one is limited in the sizes available but, since 
we are investigating in particular the topological QAH phase 
which is a translationally invariant (q = 0) instability (see 
below), we can in principle use any finite-size lattice, even 
though some of them do not have the full space-group sym¬ 
metry. On the other hand, when looking at competing charge 
instabilities, it is crucial to consider only clusters compatible 
with such orderings and we will devote some discussion on 
this topic in the next subsection. 


B. Finite lattices 

Since we plan to provide a systematic study on several clus¬ 
ters, we have implemented lattices sizes ranging from 12 to 42 
sites, with various shapes and spatial symmetries. In particu¬ 
lar, these clusters may or may not have some reflection or ro¬ 
tation symmetries (they all have inversion symmetry). More¬ 
over, some of these lattices possess the ±K points (see the 
Brillouin zone depicted in Fig. 1(b)) where the tight-binding 
dispersion exhibits two Dirac cones, leading to a 6-fold de¬ 
generacy of the half-filled free fermion ground-state (GS) so 
that correlations require some care for these clusters. We refer 
to Appendix A for further details on these finite-size clusters. 

Our large choice of clusters provides a substantial step for¬ 
ward with respect to other recent ED studies, where results 
were obtained solely on TV = 18 and TV = 24 in Ref. 10 or on 
TV = 24 and TV = 30 in Ref. 11. 

In order to illustrate the need for a systematic study of clus¬ 
ter geometries we present some data for the ground state en¬ 
ergy per site. For instance, at vanishing V\/t, as shown in 
Fig. 2(a), clusters with the K points yield the lowest energies 
for large V^/t, i.e. are compatible with the appearing order (to 
be discussed later). 



FIG. 2. (Color online) GS energy per site eo vs V 2 /t for (a) Vi s= 0 
and (b) Vi/t = 4 on various clusters. We refer to the Appendix A 
for further details on these finite-size clusters. Lines are guide to the 
eyes. 


At Vi /t as 4 and large V 2 /U the situation becomes different 
as now clusters involving M point(s) seem to have the lowest 
energies, see Fig. 2(b). This is already a clear indication of 
a phase transition between CDW (known to be stable at small 
V 2 /t and large Vi / 1) and at least two other phases that we will 
characterise later. 

Clearly, the choice of clusters is important depending on 
the parameters and the kind of instability. It was also pointed 
out recently 12 that the choice of boundary conditions might be 
crucial too. By studying the TV = 18 cluster with open bound¬ 
ary conditions (OBC), as opposed to standard periodic bound¬ 
ary conditions (PBC), the authors of Ref. 12 have found two 
level crossings for Vz/t ~ 1.6 and V 2 /t ~ 2.9 respectively 
(for Vi =0). This was taken as an indication for a stable QAH 
phase in this region. We have checked that in this case, these 
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level crossings correspond to two low-energy states having 
opposite parities with respect to inversion symmetry. How¬ 
ever, these crossings disappear when one considers the next 
largest cluster N = 32 with OBC, so that it probably corre¬ 
sponds to a short-range feature only. We generally think that 
using PBC is more appropriate to minimize finite-size effects, 
so this will be the case here. 

Last, we would like to recall some features of the semi- 
metallic (SM) phase that exists in the absence of interactions. 
Since there is a vanishing density of states at the Fermi level, 
one needs a finite strength of a short-ranged interaction to trig¬ 
ger an instability 16,17 , so that SM phase must have a finite ex¬ 
tension in the phase diagram. 

Regarding possible gap opening mechanism of the SM 
phase, Refs 18 and 19 have listed all explicit (i.e. external) 
weak-coupling perturbations which can open a gap. In the 
spinless case considered here, the three particle-hole related 
gaps are: i) the Neel-like charge density wave, which breaks 
the A-B sublattice symmetry, ii) the Kekule distortion pattern, 
which breaks translation symmetry by adopting a tripling of 
the unit-cell of modulated bond-strengths (this order parame¬ 
ter has a real and an imaginary part, thus corresponds to two 
masses), and iii) the integer quantum Hall mass 20 , induced 
by breaking the time-reversal invariance and parity symmetry 
upon adding complex Peierls phases on next-nearest neighbor 
hoppings, without enlarging the size of the unit-cell. In ad¬ 
dition to the particle-hole gaps, there is also the possibility to 
open gaps by the addition of superconducting order parame¬ 
ters 18,21,22 , we will however not address these instabilities in 
this work. 

The model Hamiltonian (1) considered here features all the 
usual symmetries. If the semi-metallic phase is gapped out 
by interactions, then the gap opening has to happen through 
the interactions by spontaneously breaking some of the sym¬ 
metries. The case i) quoted before is a well known instabil¬ 
ity, since the Neel CDW state is an obvious strong-coupling 
ground state at large V\/t. The other instabilities ii) and iii) 
currently lack a strong coupling picture, and need to be con¬ 
firmed by numerical simulations. We note however that all 
three particle-hole instabilities have been reported in mean- 
field studies. 4-6 


C. Overview of the phase diagram 

We start by drawing the global phase diagram that summa¬ 
rizes our main findings, see Fig. 3. Its main features are the 
existence of several types of charge or bond ordering for inter¬ 
mediate to large V\ and/or V 2 interactions: Neel CDW, charge 
modulation (CM), zigzag (ZZ) phase, Neel domain wall crys¬ 
tal (NDWC), and plaquette/Kekule phase (P-K) that we will 
clarify later. The large orange region (ST*) in the upper right 
part of the phase diagram features a degeneracy at the semi- 
classical level, and it is presently unclear whether and how an 
order-by-disorder mechanism will lift the degeneracy. While 
some of these phases (CM, P-K, CDW) had been predicted us¬ 
ing mean-field studies 6 and confirmed numerically in some re¬ 
gions 10,11 , the others (including NDWC and ST* phase for re¬ 


pulsive interactions and the ZZ phase for attractive V\) had not 
been advocated before. Note already that the plaquette/Kekule 
(P-K) phase only exists in some bounded region for interme¬ 
diate (Vi, Vi) values, and does not extend to strong coupling. 

There is also a large region of phase separation, mostly for 
strong attractive interactions, in agreement with the results of 
Ref. 23 for V 2 = 0. Possibly superconductivity is present in 
parts of the attractive region of the phase diagram, but we did 
not focus on this instability here. 

Last but not least, we do not have any convincing evidence 
for the stability of the QAH phase, as found in similar recent 
numerical studies 10,11 but in contradiction with another nu¬ 
merical study using open boundary conditions and entangled- 
plaquette ansatz. 12 
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FIG. 3. (Color online) Phase diagram in the (Vi /t, Vz/t) parameter 
space obtained from several exact diagonalization techniques (see 
text). Dashed lines represent the classical transition lines, see Fig. 4. 
The semi-metal, which is the ground-state for non-interacting spin¬ 
less fermions, has a finite extension in the phase diagram because of 
its vanishing density of states at the Fermi level. We will argue in the 
remainder of this article that several other phases can be stabilised for 
intermediate and/or large interactions: Neel CDW, plaquette/Kekule 
(P-K), Neel domain wall crystal (NDWC), zigzag (ZZ) phase, and 
charge modulation (CM). The region (ST*) is degenerate at the semi- 
classical level, and it is presently unclear whether and how an order- 
by-disorder mechanism will lift the degeneracy. Note also the large 
region of phase separation mostly in the attractive quadrant. Filled 
symbols correspond to numerical evidence (using level spectroscopy 
or measurements of correlations, see Sec. IV) obtained mostly on a 
N = 24 cluster which contains the most important points in its Bril- 
louin zone and features the full lattice point group symmetry of the 
honeycomb lattice. Star symbols denote likely first order transitions, 
witnessed by level crossings on the same cluster. Our numerical re¬ 
sults do not support any region of topological QAH phase. 
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We will now turn to the presentation of various numerical 
data and considerations that we have used to come up with 
this global phase diagram. 
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FIG. 4. (Color online) Classical phase diagram (a) and qualitative first order in t/(V 1 , V 2 ) phase diagram (b). The "Zigzag*" and the "Stripy*" 
phases feature a nontrivial ground state degeneracy (hence the "*" suffix in "Zigzag*" and "Stripy*"). The three points V 2 = 1,14= ±4 
and V 2 — 1, Vi = 0 feature an extensive ground state degeneracy at the classical level. Upon including the first order correction due to the 
finite hopping t, two of these points spawn new phases. The V 2 = 1, Vi = 0 point develops into the charge modulation (CM) phase, while 
the V 2 = 1, Vi = +4 point broadens into a novel "Neel domain wall crystal" (NDWC), which is sketched in Fig. 12. All regions and lines in 
(b) beyond the CM and NDWC phases have no first order (in t) quantum corrections. Note that the radial direction in the two panels has no 
physical relevance, i.e. it is not parametrizing t/V 1 , 2 . 


III. LIMIT OF LARGE INTERACTIONS (ISING LIMIT) 

The Hamiltonian (1) reduces to a free fermion problem in 
the absence of interactions, yielding a semi-metallic state with 
two Dirac points at half filling. In the opposite limit t = 0, the 
model reduces to a (generally frustrated) classical Ising model 
with competing antiferromagnetic (for Vi .2 > 0) or ferromag¬ 
netic (for Vi ,2 < 0) nearest and next-nearest neighbor inter¬ 
actions. Since the phase diagram at finite but small hopping is 
expected to be influenced by this (frustrated) Ising limit, we 
find it worthwhile to explore the classical phase diagram first. 

Our results are summarized in the phase diagrams sketched 
in Fig. 4 where we have considered the pure classical case 
t = 0 (left panel), as well as the perturbative effect of a 
small t (right panel). We will discuss these two cases and 
their phases in the following. We will exchangeably use the 
V \, V 2 notation or the 9 notation: V\ = cos 0,V 2 = sin 6, with 
0 e [0,2t r). 



FIG. 5. (Color online) Sketch of a Hamiltonian unit which renders 
the classical Vi—V 2 problem "frustration free". 


A. Classical limit: t = 0 

First we point out that this classical Hamiltonian can be 
written in a "frustration free" form, when the Hamiltonian is 
decomposed into parts acting on triangles spanned by three 
V 2 bonds and their central site (this object can alternatively 
be considered as a tripod, i.e. a site plus its three nearest V\ 
neighbors), whose three V\ couplings are counted only with 
half their strength (this is necessary because each such bond 
belongs to two different V 2 triangles, while the V 2 bonds be¬ 
long to a unique triangle), see Fig. 5. In this decomposi¬ 
tion each such tripod acts on 4 sites, and depending on 6 the 
Hamiltonian on a tripod has a different number of local ground 
states. We have numerically verified that all ground states of 
the full problem are simultaneous ground states of each tripod 
term, thus the characterization "frustration free". 

By finite cluster ground state enumeration techniques (with 
and without exploiting the frustration free property) we have 
established the classical phase diagram displayed in the left 
panel of Fig. 4. For V\ = 0, V 2 = 1 (0 = 7t/2), the 
system decouples into two interpenetrating triangular lattices 
with antiferromagnetic interactions. The ground state man¬ 
ifold is thus spanned by the product of the ground spaces 
of each triangular sublattice, and is macroscopically degen¬ 
erate (i.e. featuring an extensive ground state entropy). For 
arctan(l/4) < 9 < n /2 the ground state manifold is much 
reduced, but still grows with system size, cf. Tab. II in Ap¬ 
pendix B. On clusters featuring one or several M point(s) in 
the Brillouin zone, some of these states can be understood as 
pristine stripy ordered states (sketched in the appropriate re¬ 
gion in Fig. 4), combined with line defects, as illustrated in 
Fig. 7. It is however not clear to us whether all states in the 
ground space manifold on all lattices can be generated along 
these ideas. At 6 = arctan(l/4) we find again a single point 
with a rapidly growing ground state degeneracy. For — 7r/2 < 
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6 < arctan(l/4) we find the two Neel ground states, which 
are obvious ground states at 9 = 0 , corresponding to nearest 
neighbor interactions only. For i r — arctan(l/4) < 6 < 3it/2 
we find a region of phase separation, where the system wants 
to be either empty or completely filled with fermions. At 
0 = it — arctan(l/4) (V\ = —4, V 2 = 1) we have another 
point with an extensive ground state degeneracy, albeit with a 
smaller number than on the reflected side 6 = arctan(l/4). 
Finally the region it/2 < 9 < it — arctan(l/4) hosts a degen¬ 
erate set of states, which seem to be related to zigzag charge 
configurations (sketched in the appropriate region in Fig. 4), 
with some defect structures in addition. 

The explicit ground space counting for a number of dif¬ 
ferent clusters is listed for further reference in Tab. II in Ap¬ 
pendix B. 

B. Perturbed classical limit: t = 0 + 

After having established the classical phase diagram, we 
include the effect of a finite hopping amplitude t > 0 on the 
classical ground states. We limit ourselves here to perturbative 
arguments (mostly first order in t), and check these predictions 
by full scale ED in the following sections. 

We proceed by discussing the fate of the regions in the clas¬ 
sical phase diagram, and then analyze the behavior on the de¬ 
generate points. 

1. Neel CDW state 

We start with the least degenerate region, — it/2 < 9 < 
arctan(l/4), where there are only two Neel ground states. 
These two states are symmetry related (and thus cannot be 
split by diagonal terms) and it is not possible to connect the 
two states by a finite order in perturbation theory in the ther¬ 
modynamic limit. So based on perturbation theory this phase 
is expected to have a finite region of stability when t 7 ^ 0 . 

2 . Phase separation region 

The phase separation region it — arctan(l/4) < 0 < 3it/2 
is also stable at finite t 7 ^ 0 because both the empty and the 
full system are inert under the action of the hopping term. 

3. Zigzag state (ZZ) 

Next we consider the extended region it /2 < 6 < it — 
arctan(l/4), featuring zigzag ground states and some types 
of domain walls. To first order in t/V 1,2 it is not possible to 
hop and stay in the ground state manifold. Exact diagonaliza- 
tions of the full microscopic models on various clusters (see 
Sec. IV) with true classical ground states reveal however an 
order-by-disorder mechanism, where a diagonal term in high- 
order perturbation theory seems to lift the degeneracy and se¬ 
lect the pristine zigzag state, see Fig. 6 . This state is six-fold 


V 1 = - 8 , V 2 =6 



FIG. 6. (Color online) Splitting of the ground state degeneracy at 
Vi/t = — 8, V 2 /t = 6, yielding the two or six zigzag states as the 
ground states. Black (red) symbols denote the two different particle- 
hole symmetry sectors. All levels are singly degenerate unless the 
(spatial symmetry) degeneracy is indicated by a blue number. The 
perturbatively generated diagonal terms splits the ground space into 
symmetry related families, whereby the two or six zigzag states be¬ 
come lowest in energy. Note the very small energy difference, indi¬ 
cating a high-order perturbative effect. 


degenerate. A similar state translated to spin language appears 
in the Heisenberg-Kitaev model on the honeycomb lattice . 24 

4. Stripy* region 

Finally we consider the extended region arctan(l/4) < 
0 < 7t/2, featuring stripy ground states and some degree of 
domain walls. To first order in t/V 1,2 it is not possible to hop 
and stay in the ground state manifold. However on a finite 
cluster there is a finite order at which the moves sketched in 
Fig. 7 become possible. As the order at which this process 
first occurs diverges with the linear extent of the cluster, such 
tunneling events seem irrelevant at small t/V 1,2 in the ther¬ 
modynamic limit. As we currently do not fully understand the 
structure of the ground space manifold, we can not rule out the 
appearance of other perturbative processes connecting ground 
states without wrapping around the torus. 

In order to shed some light on this issue, we diagonalize 
the full Hamiltonian at Vi/t = 8, V^/t = 6 for a number 
of samples. The resulting spectra are shown in Fig. 8 . In 
contrast to the zigzag case studied before, here we do not ob¬ 
serve a clearly structured order-by-disorder selection of a new 
ground state. The results for these system sizes can be ra¬ 
tionalized as the spectrum of an abstract tight binding model, 
where the sites are the different ground states in the classi¬ 
cal ground space, while the hopping matrix elements are of 
the cluster wrapping type just discussed. As the perturba¬ 
tive orders of these processes strongly depend on the shape 
of the cluster, the strong correlation of the bandwidth of the 
split ground space with the formal order in perturbation the- 
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FIG. 7. (Color online) Illustration of the stripy charge ordering pattern and the line defects for the N = 24 (left panel) and N — 32 (right 
panel) clusters. The left part in each panel indicates a pristine stripy configuration, while the right part shows how another ground state can be 
reached by a one-dimensional collective transposition of particles. For the N = 24 cluster the line move wraps completely around the sample 
and connects two ground states without defect lines, while for the N = 32 the indicated shifts of particles generate a ground state with a defect 
line, where the region along the lines resembles a rotated pristine stripy state. For the N = 24 cluster this move is generated in 6th order 
perturbation theory, while the move on the N = 32 cluster is effective in 4th order. 
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FIG. 8. (Color online) Splitting of the classical ground state degen¬ 
eracy at Vi /t =5 8, V 2 /t m 6 for clusters with N = 12 (left panel) 
and N = 24, 24b, 28 and 32 sites (right panel). Black (red) symbols 
denote the two different particle-hole symmetry sectors. All levels 
are singly degenerate unless the (spatial symmetry) degeneracy is in¬ 
dicated by a blue number. The order at which processes of the type 
described in Fig. 7 happen on the clusters are indicated at the bot¬ 
tom of the frame for each system size. In contrast to the case with 
attractive Vi there is no obvious system-size independent order-by¬ 
disorder mechanism at work on the clusters considered. 


ory seems to validate our picture. It is thus presently not clear 
whether the ground state degeneracy will be lifted at a finite 
order in perturbation theory in the thermodynamic limit or 
not. The qualitative presence of a high-order diagonal term 
in the zigzag case discussed before could also carry over to 
the present stripy* case, but is masked on moderate finite size 
samples by the leading cluster-wrapping splitting terms. 


Now we turn to the points in the classical phase diagram 
which feature a massive ground state degeneracy. These are 
the points 6 = 7 t/2 (Vi = 0), 0 = arctan(l/4) (Vi = 
4, V 2 = 1) and (9 = tt - arctan(l/4) (Vi = -4, V 2 = 1) . 


5. Charge modulated (CM) state 


At Vi = 0 we have two independent triangular sublattices 
with antiferromagnetic interactions. It is known that each state 
which has two or one charges on each V 2 -triangle is a valid 
ground state. 25 When working at half filling the ground space 
of the full problem can be obtained by enumerating all ground 
states on one triangular sublattice and then summing up the 
product of the two sublattice Hilbert spaces over all particle 
bipartitions such that the total number of particles amounts to 
half filling (i.e. N/ 2). 

The effect of a finite t/V 2 is to hop particles from one trian¬ 
gular sublattice to the other one. There are indeed many states 
in the ground state manifold where such first order hoppings 
are possible while staying in the ground state manifold. The 
problem thus maps to a hopping problem on an abstract lattice 
living in the constrained Hilbert space, where the sites are the 
allowed configurations. This abstract lattice is not homoge¬ 
neous and features sites with different coordination numbers. 
It is not known to us how to solve this hopping problem an¬ 
alytically. As a starting point it is instructive to identify the 
maximally flippable states (i.e. abstract sites with the largest 
coordination number) of this problem. By explicit enumera¬ 
tion we have identified the 18 uud-ddu configurations as the 
maximally flippable states (the number is independent of TV, 
given that the sample geometry is compatible with uud-ddu 
states). These states feature a threefold degenerate charge ana¬ 
log of the magnetic uud states on one sublattice, while the 
other sublattice has three possible charge analogs of the mag¬ 
netic ddu state. Another factor of two is obtained by swapping 
the role of the two sublattices. 

The question is now whether the ground state of the full ef¬ 
fective model is localized on those maximally flippable states, 
or whether the hopping network favors delocalized wave func¬ 
tions. We have numerically solved the effective model for 
samples with TV = 18, 24, 42 and 54 sites and show the low- 
lying energy spectra in Fig. 9. We note the presence of 18 low- 
lying states (indicated by the full arrow in Fig. 9), separated 
by a gap (dashed arrow) from the higher lying states. While 
this picture is not yet very sharp for TV = 18 and TV = 24, for 
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TV = 42 and TV = 54 however the bandwidth of the 18 low- 
lying states is quite small compared to the gap to the higher 
states. We interpret this spectrum as the localization of the 
low-lying eigenstates on the 18 maximally flippable states. 
This is further corroborated by the fact that these eigenfunc¬ 
tions have their maximal amplitudes on the 18 uud-ddu states. 
We have also obtained the static charge structure factor 

S(Q) = E e_i Q ' {Rk ~ Rl \n k - 1/2)(n, - 1/2) (2) 

k,l 

in the ground state of the effective model. In Fig. 10 we dis¬ 
play this quantity for TV = 24 and TV = 42, as well as for 
an equal superposition of the 18 uud-ddu states on TV = 24 
sites. Note that we plot the data in an extended Brillouin 
zone scheme. There are obvious Bragg peaks at the momenta 
K, K' = — K stemming from the tripling of the unit cell in 
each triangular sublattice. Furthermore in the plain superposi¬ 
tion of the 18 uud-ddu states a sublattice imbalance is present, 
as the densities on the two triangular sublattices are differ¬ 
ent (1/3 versus 2/3), leading to another set of Bragg peaks 
at the T* momentum. In the actual ground state of the ef¬ 
fective model this peak is renormalized substantially. If our 
advocated picture of a localization of the ground state wave 
functions on the 18 uud-ddu states holds true, then we also 
expect the Bragg peaks at T* to survive, although with a siz¬ 
able reduction of the amplitude. In a recent paper Daghofer 
and Hohenadler 11 also investigated this issue and concluded 
an absence of the sublattice imbalance based on exact diag- 
onalizations of the original Hamiltonian (1) on V = 18 and 
TV = 24. Our results plotted in the extended zone scheme 
however reveal at least the presence of a local maximum in 
S( Q) at r*, consistent with our claim of a weak sublattice 
imbalance in the thermodynamic limit. Note that this state 
would be insulating, i.e. unrelated to the presumed metallic 
pinball liquid phase on the triangular lattice. 26 

We close the discussion of the V\ = 0 case by noting that 
this point enlarges to a finite region with respect to the addi¬ 
tion of Vi when t > 0. The reason is that the ground state of 
the hopping problem in the degenerate ground state manifold 
is able to gain energy proportional to t in the order-by-disorder 
framework. When switching on V\ , one needs a finite amount 
of Vi to destabilize the CM phase, because the stripy* region 
and the zigzag phase do not gain energy from the hopping 
term at first order in t. 


6. Neel domain wall crystal (NDWC) 



18 24 42 54 


N 

FIG. 9. (Color online) Low lying energy spectrum in units of the 
hopping t for the CM phase at Vi = 0, V 2 = 1, t/Vz 1. Data 
for system sizes N = 18, 24, 42, 54 are shown. The full arrow 
denotes the bandwidth of the lowest 18 levels. Note the collapse of 
the bandwidth for N = 42 and N = 54, while the (particle-hole) gap 
above the 18 states (dashed arrows) remains of order t, as illustrated 
in the inset. 



FIG. 10. (Color online) Charge structure factor of the effective model 
at Vi =0,V2 = l,t/V 2 1 plotted in the enlarged Brillouin zone. 
The filled symbols are data for the ground state at N = 24 while 
the empty symbols are for N = 42. The dashed symbols denote the 
structure factor in an equal superposition of the 18 uud-ddu states 
for N — 24. Note the presence of small peaks at momentum T*, 
indicating a possibly weak sublattice imbalance. 


Now let us discuss the behavior at 6 = arctan(l/4). As 
noted before, this point features a sizable ground state degen¬ 
eracy at the classical level. Since - to the best of our knowl¬ 
edge - this manifold has not been studied so far, we proceed 
by numerically diagonalizing the first order effective Hamil¬ 
tonian obtained by projecting the fermionic hopping Hamilto¬ 
nian into the classical ground state manifold. In Fig. 11 we 
present the energy per site for various clusters. It appears that 
among the studied clusters the samples with TV = 24 and 36 


sites have a particularly low energy per site (TV = 28 is close, 
but its energy per site is higher than TV = 36, so we discard 
it). It seems that this low energy correlates with a high max¬ 
imal flippability of configurations on these clusters 27 , com¬ 
pared to the other clusters. In Fig. 12 we display one such 
maximally flippable state for TV = 24. Our understanding of 
this state is that it is formed of alternating strips of the two 
Neel CDW states. Along the domain walls between the two 
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FIG. 11. (Color online) Energetics of the first order effective model 
(1st order in t/Vi j2 ) at V 1 /V 2 =4. (Top) Energy per site in units of 
t for various cluster sizes. The clusters with 24 and 36 sites have a 
particularly low energy per site. (Bottom) Single-particle charge gap 
Aip/t. This gap amounts to 4 ~ 5 t for the system sizes with a low 
energy per site. 


CDW states there are many bonds which are flippable to first 
order in t/V 1 , 2 - This state is 18 fold degenerate. There are 
3 directions in which the domain walls can run, and for fixed 
orientation, there are 3 distinct bond choices for the domain 
walls. The remaining factor of two is due to the assignment 
of the two Neel CDW states. We have further checked a sam¬ 
ple with TV = 72 and again found 18 maximally flippable 
states, in agreement with our analysis. In Fig. 13 we display 
the low lying eigenstates of the effective first-order Hamilto¬ 
nian in the degenerate subspace for system sizes N = 24,36 
and 72 sites. The collapse of the lowest 18 states is clearly 
visible, providing strong evidence for symmetry breaking of 
the NDWC type. The results of the effective model simula¬ 
tions for the charge structure factor are shown in Fig. 14 and 
further corroborate this picture. 

We also expect the point V 1 /V 2 = 4 to enlarge to an actual 
phase at small but finite t/V i ? 2 , because the two neighbor¬ 
ing phases in the Ising limit (the stripy* region and the Neel 
CDW) both do not gain energy to first order int/Vi^- 


7. Vi = -4, V 2 = 1 

The ground state configurations at V\ = —4, V 2 = 1 do 
not seem to be flippable to first order in t/V \, V 2 , and we have 
therefore refrained from a detailed study of this line. 

To conclude this section we show the resulting semiclassi- 
cal qualitative phase diagram in the right part of Fig. 4. 


IV. NUMERICAL RESULTS 


In this section, we will present our systematic numerical 
study using ED and typical observables that were used to com¬ 
pute our global phase diagram in Fig. 3. We will first present 



Neel B 


Neel A 


Neel B 


FIG. 12. (Color online) The Neel Domain Wall Crystal (NDWC): 
sketch of a classical ground state at Vi =4, V 2 = 1 on the N = 
24 sample, which is maximally flippable within the classical ground 
state manifold with respect to the hopping t. The shaded regions 
denote the two Neel domains, and the orange circled bonds along the 
domain walls are flippable to first order in t. The green box indicates 
a twelve-site unit cell. 



FIG. 13. (Color online) Low lying energy spectrum in units of the 
hopping t for Vi = 4, V 2 = 1, t/U 2 <C 1. Data for system sizes 
N — 24, 36, 72 are shown. The full arrow denotes the bandwidth 
of the lowest 18 levels (N = 36 does not have rotational invariance, 
in that case we only expect 6 possible ground states). Note the col¬ 
lapse of the bandwidth for larger clusters, while the (particle-hole) 
gap above the 18 states (dashed arrows) remains approximately con¬ 
stant, as illustrated in the inset. 


some useful tools to detect phase transitions as well as to char¬ 
acterize phases. Then we will provide numerical evidence for 
the finite extension of the phases that were discussed in the 
limit of large interactions in Sec. III. Last but not least, we will 
discuss the stability of phases (with a stronger quantum char¬ 
acter) at intermediate coupling, namely the semi-metal (SM), 
the plaquette/Kekule phase and the absence of the conjectured 
topological Mott insulator phase. 
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FIG. 14. (Color online) Charge structure factor of the effective model 
at Vi =4,1/2 = 1, t/V 2 1 plotted in the enlarged Brillouin zone. 
The filled (empty) symbols are data for the ground state at N = 24 
(N = 72). The Bragg peaks appear at the X * point, c.f. Fig. 1(b). 


A. Numerical tools 

1. Energetics 

Let us start by showing some full low-energy spectra, where 
we label each eigenstate according to its quantum numbers. 28 
In Fig. 15-16, we present these spectra for a fixed V\ = 0 and 
V\/t = 4 respectively on the N = 24 cluster which has the 
full point-group symmetry Cq v and also possesses the relevant 
K and M points in its Brillouin zone (see Appendix A). 

When Vi/t = 0 and V 2 /t is small, we expect a finite region 
of SM phase with low-energy excitations at the Dirac point 
(±K). Increasing V 2 /t, there is a sudden level crossing when 
V 2 ft — 2.2. Beyond this value, if we count the lowest energy 
levels well separated from the others, we get 18 states, which 
may indicate some ordering of some kind in the thermody¬ 
namic limit (see discussion in Sec. IIIB). 

When Vi/t = 4, starting at small V 2 /t, there are two 
almost-degenerate states at the T point with opposite particle- 
hole symmetries, compatible with the Neel CDW state. In¬ 
creasing V 2 /t leads to some intermediate phase 1 < V 2 /t < 2 
with low-energy states at momentum K compatible with a 
Kekule-like (or plaquette) pattern, and finally for large V 2 /t > 
2 the emergence of 6 low-energy states with momenta T and 
M. 


B. Stability of the classical Ising-like phases 

Let us now move to a more precise characterization of the 
possible patterns occurring for large interactions. Quite gener- 
ically, large density interactions are expected to lead to some 
kind of charge ordering. In the equivalent spin language, these 
terms are of the Ising type and the resulting Ising model would 



FIG. 15. (Color online) Full spectrum vs V 2 /t at fixed Vi/t = 0 
labeled with symmetry sectors on N = 24 cluster. More precisely, 
states with momentum T can be labeled with Cq v standard point- 
group notations; with momentum X only reflection can be used to 
detect even (e) vs odd (o) states; at the K point, we use C^ v point- 
group notations; at the M point, we can use two reflections to la¬ 
bel states which are even/even (el), even/odd (e2), odd/even (ol) or 
odd/odd (o2). Open (filled) symbols denote respectively states which 
are even (odd) with respect to particle-hole symmetry. 



FIG. 16. (Color online) Same as Fig. 15 at fixed V\/t — 4 on A = 
24 cluster. For intermediate values V^/t — 1.5, low-energy states 
(indicated with an ellipse) are compatible with a plaquette or Kekule 
phase. 


correspondingly exhibit some spin ordering. For model (1), 
they were discussed in Sec. Ill and are denoted under the 
names: CDW, NDWC, stripy* region, CM or zigzag. 


1. CDW phase 

When V 2 = 0, since the honeycomb lattice is bipartite, a 
simple charge-density wave (CDW) where particles are local¬ 
ized on one sublattice is expected for V\ > 0, so that de- 
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FIG. 17. (Color online) Connected density correlation functions be¬ 
tween a reference site (open circle) and other sites for Vi = 1 and 
Vi = 2 (at fixed V 2 = 0) on N = 32 cluster. Periodic boundary 
conditions are used. Blue and red correspond to positive/negative 
correlations. 

generacy is two. Since this phase does not break translation 
symmetry but only inversion, symmetry analysis predicts that 
on a finite cluster there should be a collapse of the lowest YA 2 
and TD 1 states. 29 This is indeed what is observed in Fig. 16 
at small V^/t for instance. When V\ > 2, our numerical data 
are compatible with a vanishing gap at zero momentum (data 
not shown) because of the inversion symmetry breaking in the 
CDW phase. Note however that the system remains an insu¬ 
lator. 

A clearer picture of this CDW state is provided by density- 
density correlations. On Fig. 17, we plot these connected cor¬ 
relations obtained with a N — 32 cluster showing a clear tran¬ 
sition from a semi-metallic phase to the Neel CDW for large 
V\jt > 2 (at V 2 =0). About the critical value, let us mention 
that a recent unbiased QMC approach 14 has concluded that 
(y x !t) c = i.36. 

Based on our previous semiclassical analysis, we expect 
that the CDW phase cannot extend beyond the line V\ = 4 V 2 , 
at least towards strong coupling. 

2. Charge modulated phase 

As we have discussed in Sec. Ill, the large V 2 /t > 0 limit is 
less trivial since the classical Ising model possesses an ex¬ 
tensive ground-state degeneracy. However, one expects on 
general grounds that quantum fluctuations will provide some 
selection mechanism within these low-energy states. Our dis¬ 
cussion in Sec. IIIB 5 has led us to conclude in favor of some 
charge-modulation (CM) state which is 18-fold degenerate 
and exhibits charge imbalance between the two sublattices, 
in agreement with Refs. 6 and 10 and slightly different from 
the similar CM state without charge imbalance 11 . 

Our ground-state energy plot in Fig. 2(a) has indeed shown 
that the emergent phase is compatible with clusters having the 
K point, as is the case for the proposed tripled unit cell. More¬ 
over, the low-energy spectra on these clusters, such as N = 24 
in Fig. 15 (or similarly N = 36, data not shown), do confirm 
the presence of 18 low-energy states at large V 2 in agreement 
with the expected degeneracy. The question whether all those 


states will collapse or if there will remain a finite gap (and 
thus a smaller degeneracy) is left for future studies, but our 
effective model diagonalizations in the strong coupling limit 
indicate a complete collapse of the 18 levels, c.f. Fig. 9. 


3. Stripy* region 

For large V 2 ~ V\ > 0, the lowest energies are found on 
clusters having the M point (see Fig. 2(b)), which indicates a 
different instability. Moreover, our symmetry resolved spec¬ 
tra (see Fig. 16 for instance) predict a 6-fold degeneracy for 
large V 2 jt (with 3 states at momentum T and 3 at equivalent 
momenta M) in disagreement with both the CM phase or the 
Kekule one (see later). As discussed in Sec. IIIB 4 based also 
on ED, the fate of the stripy* region at finite t/V 1 , V 2 is not 
clear. Results for N = 32 close to the phase transitions to the 
neighboring phases still do not shed light on a possible selec¬ 
tion process among the classically degenerate ground states. 


4. Neel domain wall crystal 

Since the ratio V 1 /V 2 = 4 is a special point in the classical 
phase diagram we investigate here the phase diagram along 
the line V 2 /U with the ratio V 1 /V 2 = 4 fixed. In Fig. 18 we 
display the low energy spectrum for the N = 24 site sam¬ 
ple along this line. At very small V 2 /t we expect the sta¬ 
ble semimetallic phase emerging out of V 2 /t = 0. We then 
detect an intermediate Kekule or plaquette phase for values 
1 < V 2 /t < 2. These values are confirmed by directly com¬ 
puting kinetic correlations and checking if they match the ex¬ 
pected pattern for a plaquette phase. However for larger val¬ 
ues V 2 ft > 3 the system enters the realm of the NDWC phase, 
involving the physics discussed in subsection IIIB 6. The oc¬ 
currence of two very distinct phases beyond the semimetallic 
phase is also visible in the ground state correlations shown in 
Fig. 19. 


5 . Zizag phase 

As already discussed in Sec. IIIB, ED with large negative 
Vi/t and positive V 2 /t provide spectral evidence for a quan¬ 
tum order-by-disorder selection of a pristine zigzag state, as 
shown for instance in Fig. 6. For the clusters having the 3 
M points in their Brillouin zone (see Appendix A), we do 
observe 6-fold quasi degenerate ground-states. By comput¬ 
ing connected density correlations on any of these states, we 
obtain a pattern compatible with zigzag state. Indeed, a sim¬ 
ple direct computation of connected density correlations on 
the zigzag state leads to values ±1/12 or ±1/4 depending 
on distances between sites, which agrees to a very high accu¬ 
racy with the exact results computed on N = 24 or N = 32 
ground-states in this region, see Fig. 20. 

Having clarified the phases at large interaction strengths, 
which are connected to the classical limit, we now con- 
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FIG. 20. (Color online) Connected density correlations computed on 
N = 24 (left panel) and N = 32 (right panel) cluster for Vi /t = 
—8 and Vh/t — 6. Blue and red correspond to positive/negative 
correlations. 


FIG. 18. (Color online) Same as Fig. 15 as a function of V 2 /£ at fixed 
V 1 /V 2 = 4 on N = 24 cluster. The orange boxes denote the two 
parameter values for which the correlations shown in Fig. 19 were 
calculated. 


V 1 /t = 4 
V 2 /t = 1 



V x jt = 40 
V 2 /t = 10 



kinetic energy 
correlations 


• .« 



Kekule/Plaquette 

phase 



density 

correlations 


Neel Domain 
Wall Crystal 


FIG. 19. (Color online) Kinetic energy and density correlations 
computed on N = 24 cluster along the V 1 /V 2 — 4 line, i.e. for 
(Vi /t, V 2 /t) respectively equal to (4,1) and (40,10). 


sider the intermediate interaction regime where new quantum 
phases are expected to emerge. 


C. Emergence of a resonating quantum plaquette phase or 
Kekule one 

In Refs. 6 and 10, it has been argued using mean-field or 
N = 18 ED respectively that a Kekule pattern can be stabi¬ 
lized when V\ and V 2 are close. This would correspond to a 
3-fold degeneracy 30 31 of the ground-state (more precisely one 
state at T and two states at ±K). 

However, while the Kekule phase was seen to remain stable 
up to large V\ ~ V 2 in Refs. 6 and 10, our numerical data 


on large clusters do not confirm these findings (see Sec. IV B) 
since it is replaced either by the stripy* region or the NDWC 
phase. However, we do confirm its existence, but for interme¬ 
diate values only. 

A first evidence is given by the low-energy spectrum of 
N = 24 at fixed V\/t = 4 (see Fig. 16) where in some inter¬ 
mediate V 2 /t range, the lowest energies are compatible with 
this plaquette/Kekule phase. A second one is provided via 
a direct computation of the (connected) kinetic energy corre¬ 
lations on a larger N = 42 cluster in Fig. 21, which are in 
perfect agreement with the expected pattern. 30 



FIG. 21. (Color online) Kinetic energy correlations at fixed Vi/t — 8 
and V 2 /1 — 2 on N — 42 cluster. Periodic boundary conditions 
are used. Blue and red correspond to positive/negative correlations. 
Width is proportional to the correlation. 

There is a small caveat here, since it is rather difficult to dis¬ 
tinguish a static Kekule pattern from a plaquette crystal, where 
the three strong bonds on a hexagon resonate. Indeed both 
phases break the same symmetries and have the same degen¬ 
eracy so that a clear identification is often difficult. 30,32,33 Fet 
us mention that for S = 1/2 Heisenberg model on the hon- 
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eycomb lattice with competing Ji, J 2 (and J3) interactions, 
there exists an intermediate region with a clear plaquette char¬ 
acter . 30,34-36 Moreover, in this case, the plaquette phase seems 
connected to the ordered phase through a possibly continu¬ 
ous phase transition . 30,35 In our spinless fermionic model, we 
have not fully clarified the nature of this intermediate phase, 
but this would be an interesting project as well as the nature of 
the transition both to the semi-metallic and the Neel phases, 
which is unfortunately beyond what can be done using ED 
technique, but which might be fermionic examples of phase 
transitions beyond the Landau paradigm. 


D. Stability of the topological phase 

The proposal made by Raghu et al 4 is that, when V 2 > 0 
dominates, the system undergoes a transition to a topolog¬ 
ically distinct phase that spontaneously break time-reversal 
symmetry. A simple picture of this phase can be obtained by 
drawing circulating currents on a given sublattice. In order to 
check this picture, we have computed current-current corre¬ 
lations on a given sublattice. Note that these correlations are 
strictly zero for all distances when V\ = V 2 = 0. They are 
plotted for different V 2 in Fig. 22 in the case of N = 24 lat¬ 
tice and at fixed V\ — 0 where the instability is expected to be 
the strongest. 37 We have used the expected q = 0 orientation 
pattern of the QAH phase, i.e. where all currents would be 
in the clockwise direction on each hexagon for instance. We 
see on Fig. 22 that when we increase V^/t, the signal starts by 
increasing in amplitude, and then some defects (in red) start to 
appear on the cluster, in particular in the case where QAH is 
supposed to be the strongest around V 2 /t ~ 2. Let us empha¬ 
size again that QAH is perfectly compatible with all clusters 
(since it is a q = 0 instability) and as a result no frustration 
is expected. Therefore, the appearance of defects reveal some 
competing phases. Still the overall pattern is almost perfect, 
which is remarkable and does not occur in the V\ only case 
for instance. For even larger V 2 , the signal suddenly disap¬ 
pears (see Figs. 23 and 24 for several clusters, other data not 
shown), signalling the entrance in the CM phase (see above). 
Given that the semimetallic phase should have a finite exten¬ 
sion for small V^/t, this constrains the putative QAH phase in 
a small region around — 2. 

In order to estimate the finite-size effects, we have per¬ 
formed extensive simulations on a large number of samples 
(data not shown). We can clearly identify some clusters on 
which the patterns are far from being perfect (i.e. have many 
defects), which can be related either to lower symmetry, or to 
poorly two-dimensional aspect ratio. Indeed, when looking at 
clusters having Cq symmetry at least, we do observe nice cur¬ 
rent correlations with no or few defects at small V^/t. There¬ 
fore, there are definitely correlations compatible with QAH 
phase, but one needs to investigate whether these are short- 
range features only or not. 

In order to do so, we have to measure some order parame¬ 
ter. We have chosen the current structure factor, i.e. the sum 
of current correlations with the signs taken from the QAH pat¬ 


tern: 



«ii» 


where the sum runs over all A /5 bonds on a given sublattice 38 
and the £{j = ±1 correspond to the expected QAH orienta¬ 
tion (all currents clockwise or anticlockwise depending on the 
choice of reference bond (iojo))- When using clusters with 
less symmetry, we average over the inequivalent directions of 
the reference bond. This is indeed a valid order parameter 
and long-range order requires that it diverges in the thermo¬ 
dynamic limit as ~ N. Overall, we see in Fig. 25 a clear 
increase of current correlations for moderate V 2 , indicating 
that such an interaction indeed favours current formations, in 
agreement with mean-field prediction. There are of course 
non-trivial finite-size effects due to the fact that not all clus¬ 
ters have the same symmetry (see Appendix A). 

In the inset of Fig. 25, we plot this structure factor divided 
by the number of terms for 1 < V 2 /t < 2 where the QAH 
seems to be the strongest from our data, and also from the 
mean-field estimate of the QAH region 4 : 1.5 ;$ V 2 2.5. 
For clusters without the K point, the behaviour of this cur¬ 
rent structure factor has some non-monotonic size dependence 
due to the fact that clusters have different shapes, but its ab¬ 
solute value is quite small, and finite-size extrapolations are 
compatible with a vanishing value. We also plot separately 
the data points for clusters having the K point (which have a 
larger signal due to the degenerate ground state at half filling 
at V 2 = 0) but if we focus on N = 24 and N = 42 which are 
the best two-dimensional clusters of this kind (good aspect ra¬ 
tio and Cq symmetry at least), data are also compatible with 
the absence of long-range order in the thermodynamic limit. 
Our results are in agreement with other recent numerical stud¬ 
ies 10,11 but in contrast with another work . 12 We hope that fu¬ 
ture progress in simulations of two-dimensional strongly cor¬ 
related systems will allow to confirm the absence of a QAH 
phase for this microscopic model, but will be able to stabilise 
it in a larger parameter space. Indeed, it is obvious that the 
semi-metallic state has a strong response in the QAH channel, 
so that if other instabilities could be prevented (using addi¬ 
tional interactions to frustrate them), there remains an inter¬ 
esting possibility to stabilise QAH phase in a truly long-range 
ordered fashion. 


V. CONCLUSION 

We have investigated a model of spinless fermions on hon¬ 
eycomb lattice with (generally competing) density-density in¬ 
teractions using extensive numerical exact diagonalizations on 
various clusters. Based on a careful analysis, our main result 
is summarized in the global phase diagram in Fig. 3. On the 
one hand, we have clarified the phase diagram in the strong 
coupling regime (\Vi^\ t) using some non-trivial results 
on the classical model and some quantum order-by-disorder 
arguments. On the other hand, we have used numerical simu¬ 
lations of the full quantum model on finite clusters. As a re¬ 
sult, for intermediate or large interactions, we have provided 
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FIG. 22. (Color online) Current-current correlations on a given sublattice between a reference bond (black) and other bonds for V 2 /t = 1, 
2 and 3 (from left to right) with fixed Vi = 0 on N = 24 sample. Periodic boundary conditions are used. Blue and red correspond to 
positive/negative correlations according to the orientations expected in the QAH phase (taken here to be clockwise on all hexagons). Width is 
proportional to the correlation. 



FIG. 23. (Color online) Same as Fig. 22 for V± = 0 and V^/t — 1 on various clusters having Cq symmetry. From left to right: N = 26, 32, 
38 and 42. 





FIG. 24. (Color online) Same as Fig. 23 for Vi = 0 and Vz/t — 2. 


evidence for several crystalline phases: Neel CDW, CM, pla- 
quette/Kekule, Neel domain wall crystal, the stripy* region 
and the zigzag phase, some of which had not been considered 
yet. 


When considering the role of V 2 >0 interaction alone, it 
had been suggested in the literature 4 that it could stabilize 
a topological phase, which has triggered an intense activity. 
We have measured the corresponding (charge) current corre- 
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FIG. 25. (Color online) Current structure factor J q as a function of 
V 2 It for several lattices. Inset: scaling of this structure factor divided 
by the number of bonds Nb = 3N/2 — 11 vs 1/N for several V 2 /t. 
Open symbols are used for clusters without the K point. Hatched 
symbols are used for the clusters N — 24 and N — 42 which have 
the K point as well as Cq symmetry at least. All data points are 
compatible with a vanishing signal in the thermodynamic limit. 

lations, and while there is definitely some increase due to the 
V 2 interaction, our finite-size analysis indicate that the topo¬ 
logical phase is not present in the thermodynamic limit, and 
there is a direct transition between the semi-metallic phase and 
a charge-modulated one. However, it is very interesting to no¬ 
tice that V 2 interaction is able to create short-range current 
correlations, so that some additional ingredients could possi¬ 
bly stabilize it. Future work along these lines should also con¬ 
sider the spinful case where a topological spontaneous quan¬ 
tum spin hall phase with spin currents has been proposed. 4 

As a last remark, let us remind that other models are known 
to exhibit Dirac cones and thus could potentially also host 
such a topological phase: for instance, on the square lattice 
7r-flux model at half-filling, a topological insulating phase 
is predicted using mean-field technique 5 but numerical study 
have not confirmed its stability. 39 It is also noteworthy that a 
generalized tight-binding model on the kagome lattice always 
contain Dirac cones at some particular filling 40 , which leaves 
room to investigate the exciting possibility to stabilize a topo¬ 
logical phase with interactions. 

Note added: While finalizing this manuscript we became 
aware of concurrent iDMRG work 41 with consistent results 
regarding the presence and absence of the various phases in 
the repulsive regime. 


port. AML was supported by the FWF (I-1310-N27/DFG 
FORI 807). 

Appendix A: Lattice geometries 


Since we are using several kinds of lattices, we give their 
definitions in Table I using a and b translations to define the 
torus with periodic boundary conditions. 


Name 

a 

b 

sym. group 

K (2) 

M (3) 

X(6) 

12 

(3,3^3) 

(3,— s/3) 

z 2 

no 

yes(l) 

no 

18 

(6,0) 

(3,3a/3) 

Cq v 

yes 

no 

no 

24 

(6,2 a/3) 

(0, 4a/3) 

Cq v 

yes 

yes 

yes 

24b 

(4,4^3) 

(-4,2^3) 

z 2 

no 

yes(l) 

no 

26 

(7, V3) 

(2,4^3) 

c 6 

no 

no 

no 

28 

(7, V3) 

(0,4a/3) 

Z2 

no 

yes(1) 

no 

30a 

(3, 5a/3) 

(-3,5a/3) 

C 2 X C2 

yes 

no 

no 

30b 

(5,3^3) 

(-5,3^3) 

C 2 x C 2 

no 

no 

no 

32 

(8,0) 

(4,4^3) 

Cq v 

no 

yes 

no 

34 

(9,a/3) 

(2.4V3) 

Z 2 

no 

no 

no 

36 

(6,0) 

(6,6a/3) 

C 2 X C 2 

yes 

yes (1) 

yes (2) 

38 

(8,2 a/3) 

(1, 5a/3) 

C 6 

no 

no 

no 

42 

(9, a/3) 

(3,5 a/3) 

C 6 

yes 

no 

no 

54 

(9,3a/3) 

(0,6a/3) 

Cq v 

yes 

no 

no 

72 

(8,0) 

(8,8a/3) 

Cq v 

yes 

yes 

yes 


TABLE I. Finite lattices studied in this work. Listed are the number 
of spins N ; the basis vectors a, b in the plane; the symmetry point 
group; presence of the K point; M point; X point. 


Appendix B: Counting of Ising Groundstates 

In this section we provide a table with the dimensions of the 
classical ground spaces at half filling for the clusters described 
in Appendix A. 
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N s 

- 4/5 

-4/5<Vi<0 

0 

0<Vi<4/5 

4/5 

12 

- 

8 

- 

8 

- 

18 

- 

- 

666 

0 

38 

24 

96 

6 

5’526 

6 

446 

24b 

- 

12 

- 

12 

- 

26 

- 

- 

3’042 

0 

340 

28 

- 

- 

8’964 

2 

858 

30a 

- 

- 

30’618 

0 

692 

30b 

- 

- 

11’300 

0 

888 

32 

630 

42 

1*764 

42 

1*520 

36 

- 

- 

230’014 

14 

3’454 

38 

- 

- 

- 

- 

5’208 

42 

- 

- 

1’211’004 

- 

8’880 

54 

- 

- 

59’711’598 

- 

- 

72 

- 

186 

- 

186 

8’705’390 


TABLE II. Ground state degeneracy in the classical limit t — 0 at 
half filling. 0 denotes a ground state energy on this cluster that is 
higher than on other clusters. - indicates that the ground space has 
not been explored for this instance. The column labels Vi and we use 
the parametrization V 2 = 1 — Vi for Vi > 0 and V2 = 1 — | Vi | for 
Vi < 0. 
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